Log-Laplace Distribution — a log-symmetric, heavy-tailed model on \((0,\infty)\)

The log-Laplace distribution (SciPy: scipy.stats.loglaplace) is a positive distribution whose logarithm is Laplace (double-exponential).

It’s a useful choice when you want a distribution on positive magnitudes with:

• symmetry on a multiplicative (log) scale,

• power-law tails (heavier than lognormal), and

• a simple inverse-CDF sampler.

What you’ll learn

• How to define the log-Laplace PDF/CDF and connect it to the Laplace distribution via a change of variables.

• Which moments exist (and why some don’t), including closed forms for mean/variance/skewness/kurtosis when they do.

• A NumPy-only inverse-CDF sampler and how to validate it.

• Simple inference for the shape parameter: MLE, exact confidence intervals, and a conjugate Bayesian update.

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import stats

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")  # CKC convention

np.set_printoptions(precision=5, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & Classification

• Name: Log-Laplace distribution (loglaplace)

• Type: Continuous

• Support (standard): \(x \in (0,\infty)\)

• Parameter space (standard): shape parameter \(c>0\)

SciPy uses an additional location/scale transform:

\[X = \text{loc} + \text{scale}\,Y, \quad Y \sim \mathrm{LogLaplace}(c),\]

where scale > 0 and the support becomes \(x > \text{loc}\).

We write (standard form):

\[X \sim \mathrm{LogLaplace}(c).\]

2) Intuition & Motivation

2.1 What it models

A log-Laplace random variable can be written as

\[X = \exp(Z), \quad Z \sim \mathrm{Laplace}(0, 1/c).\]

So additive Laplace noise in log-space becomes multiplicative, heavy-tailed noise on the original scale.

A key symmetry is:

\[\log X \text{ is symmetric about } 0 \quad\Longleftrightarrow\quad X \text{ is symmetric about } 1 \text{ on the multiplicative scale.}\]

In fact, \(X\) and \(1/X\) have the same distribution.

2.2 Typical real-world use cases

• Positive data with occasional extreme values: file sizes, response times, claim sizes.

• Multiplicative error models: \(Y = \theta \cdot X\) where \(\log X\) has Laplace-like “spiky + outlier-prone” behavior.

• Robust alternatives in log-space: if \(\log Y\) is better modeled by Laplace than Normal (heavier tails), then \(Y\) may be log-Laplace rather than lognormal.

2.3 Relations to other distributions

• Laplace: \(\log X \sim \mathrm{Laplace}(0, 1/c)\).

• Lognormal: \(\log X\) Normal vs Laplace (log-Laplace has sharper peak and heavier tails).

• Pareto-like tails: for \(x\to\infty\), \(f(x)\propto x^{-(c+1)}\).

• Invariance: \(X \stackrel{d}{=} 1/X\) (log-symmetry).

3) Formal Definition

3.1 PDF

The standard log-Laplace density can be written in a compact way:

\[f(x; c) = \frac{c}{2x}\,\exp\bigl(-c\,|\log x|\bigr), \quad x>0,\; c>0.\]

Equivalently, as a piecewise power law (useful for intuition):

\[\begin{split} f(x;c) = \begin{cases} \frac{c}{2}\,x^{c-1}, & 0 < x < 1,\\[4pt] \frac{c}{2}\,x^{-c-1}, & x \ge 1. \end{cases} \end{split}\]

3.2 CDF

\[\begin{split} F(x;c) = \begin{cases} 0, & x \le 0,\\[4pt] \frac{1}{2}x^{c}, & 0 < x < 1,\\[4pt] 1 - \frac{1}{2}x^{-c}, & x \ge 1. \end{cases} \end{split}\]

3.3 Quantile function (inverse CDF)

This is especially convenient for sampling:

\[\begin{split} F^{-1}(p;c)= \begin{cases} (2p)^{1/c}, & 0<p<\tfrac12,\\[4pt] \left(\frac{1}{2(1-p)}\right)^{1/c}, & \tfrac12\le p<1. \end{cases} \end{split}\]

def _validate_c(c: float) -> float:
    c = float(c)
    if c <= 0:
        raise ValueError("c must be > 0")
    return c


def loglaplace_pdf_std(x: np.ndarray, c: float) -> np.ndarray:
    """Standard log-Laplace PDF on (0, inf) with shape c>0."""
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)

    out = np.zeros_like(x, dtype=float)
    mask = x > 0
    xm = x[mask]

    # Stable expression: f(x)= (c/(2x))*exp(-c*|log x|)
    out[mask] = 0.5 * c / xm * np.exp(-c * np.abs(np.log(xm)))
    return out


def loglaplace_logpdf_std(x: np.ndarray, c: float) -> np.ndarray:
    """Standard log-Laplace log-PDF."""
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)

    out = np.full_like(x, -np.inf, dtype=float)
    mask = x > 0
    xm = x[mask]

    out[mask] = np.log(c) - np.log(2.0) - np.log(xm) - c * np.abs(np.log(xm))
    return out


def loglaplace_cdf_std(x: np.ndarray, c: float) -> np.ndarray:
    """Standard log-Laplace CDF."""
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)

    out = np.zeros_like(x, dtype=float)
    mask = x > 0
    xm = x[mask]

    out[mask] = np.where(
        xm < 1.0,
        0.5 * np.power(xm, c),
        1.0 - 0.5 * np.power(xm, -c),
    )
    return out


def loglaplace_ppf_std(p: np.ndarray, c: float) -> np.ndarray:
    """Standard log-Laplace inverse CDF."""
    p = np.asarray(p, dtype=float)
    c = _validate_c(c)
    if np.any((p <= 0) | (p >= 1)):
        raise ValueError("p must be in (0, 1)")

    return np.where(
        p < 0.5,
        np.power(2.0 * p, 1.0 / c),
        np.power(1.0 / (2.0 * (1.0 - p)), 1.0 / c),
    )


def loglaplace_rvs_numpy(
    c: float, size: int | tuple[int, ...], rng: np.random.Generator
) -> np.ndarray:
    """NumPy-only sampling via inverse CDF (standard form)."""
    c = _validate_c(c)

    u = rng.random(size)
    # Keep u in (0,1) to avoid returning exactly 0 or inf from floating endpoints.
    u = np.clip(u, np.finfo(float).tiny, 1.0 - np.finfo(float).eps)

    return np.where(
        u < 0.5,
        np.power(2.0 * u, 1.0 / c),
        np.power(1.0 / (2.0 * (1.0 - u)), 1.0 / c),
    )


def loglaplace_pdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    """Location/scale version: X = loc + scale*Y, Y ~ LogLaplace(c)."""
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)
    scale = float(scale)
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale
    return loglaplace_pdf_std(z, c) / scale


def loglaplace_logpdf(
    x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0
) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)
    scale = float(scale)
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale
    return loglaplace_logpdf_std(z, c) - np.log(scale)


def loglaplace_cdf(x: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)
    scale = float(scale)
    if scale <= 0:
        raise ValueError("scale must be > 0")

    z = (x - loc) / scale
    return loglaplace_cdf_std(z, c)


def loglaplace_ppf(p: np.ndarray, c: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    c = _validate_c(c)
    scale = float(scale)
    if scale <= 0:
        raise ValueError("scale must be > 0")

    return loc + scale * loglaplace_ppf_std(p, c)


# Quick cross-check against SciPy (standard form: loc=0, scale=1)
x_test = np.array([1e-3, 0.2, 0.9, 1.0, 2.0, 10.0])
c_test = 2.5
dist = stats.loglaplace(c_test)

print("max |pdf - scipy|:", np.max(np.abs(loglaplace_pdf_std(x_test, c_test) - dist.pdf(x_test))))
print("max |cdf - scipy|:", np.max(np.abs(loglaplace_cdf_std(x_test, c_test) - dist.cdf(x_test))))

max |pdf - scipy|: 3.7947076036992655e-19
max |cdf - scipy|: 0.0

4) Moments & Properties

A convenient fact is that the log-Laplace has power-law tails, so not all moments exist.

4.1 Raw moments

For real \(k\) with \(-c < k < c\),

\[\mathbb{E}[X^k] = \frac{c^2}{c^2 - k^2}.\]

In particular:

• Mean exists iff \(c>1\): $\(\mathbb{E}[X] = \frac{c^2}{c^2 - 1}.\)$

• Second moment exists iff \(c>2\): $\(\mathbb{E}[X^2] = \frac{c^2}{c^2 - 4}.\)$

4.2 Variance

For \(c>2\):

\[\mathrm{Var}(X) = \frac{c^2(2c^2+1)}{(c^2-4)(c^2-1)^2}.\]

4.3 Skewness and kurtosis

These require higher moments:

• Skewness exists for \(c>3\):

\[ \gamma_1 = \frac{2(15c^4+7c^2+2)\,\sqrt{c^2-4}}{c\,(c^2-9)\,(2c^2+1)^{3/2}}. \]

• Excess kurtosis exists for \(c>4\):

\[ \gamma_2 = \frac{6\bigl(2c^{10}+138c^8-615c^6-449c^4-132c^2-24\bigr)}{c^2\,(c^2-16)\,(c^2-9)\,(2c^2+1)^2}. \]

4.4 MGF / characteristic function

Because the right tail behaves like \(x^{-(c+1)}\), the moment generating function

\[M_X(t)=\mathbb{E}[e^{tX}]\]

diverges for any \(t>0\) (too much mass in the far right tail for an exponential weight).

For \(t<0\), the Laplace transform exists and can be written using incomplete gamma functions:

\[ \mathbb{E}[e^{tX}] = \frac{c}{2}\Bigl[(-t)^{-c}\,\gamma(c,-t) + (-t)^{c}\,\Gamma(-c,-t)\Bigr], \quad t<0. \]

The characteristic function \(\varphi_X(\omega)=\mathbb{E}[e^{i\omega X}]\) exists for all real \(\omega\) (bounded integrand), and admits an analogous representation with complex arguments.

4.5 Entropy

Using \(\log X \sim \mathrm{Laplace}(0,1/c)\) and the entropy change-of-variables rule,

\[h(X) = 1 + \log\Bigl(\frac{2}{c}\Bigr).\]

(Here \(h\) is differential entropy.)

4.6 Other handy properties

• Median: \(\mathrm{median}(X)=1\).

• Log-moments: \(\mathbb{E}[\log X]=0\) and \(\mathrm{Var}(\log X)=2/c^2\).

• Reciprocal invariance: \(X \stackrel{d}{=} 1/X\).

def loglaplace_raw_moment(k: float, c: float) -> float:
    """Return E[X^k] for the standard log-Laplace, when it exists (|k|<c)."""
    c = _validate_c(c)
    k = float(k)
    if not (-c < k < c):
        return np.inf
    return (c * c) / (c * c - k * k)


def loglaplace_mean(c: float) -> float:
    c = _validate_c(c)
    return loglaplace_raw_moment(1.0, c) if c > 1 else np.nan


def loglaplace_variance(c: float) -> float:
    c = _validate_c(c)
    if c <= 2:
        return np.nan
    c2 = c * c
    return c2 * (2 * c2 + 1) / ((c2 - 4) * (c2 - 1) ** 2)


def loglaplace_skewness(c: float) -> float:
    c = _validate_c(c)
    if c <= 3:
        return np.nan
    c2 = c * c
    num = 2 * (15 * c2 * c2 + 7 * c2 + 2) * np.sqrt(c2 - 4)
    den = c * (c2 - 9) * (2 * c2 + 1) ** 1.5
    return num / den


def loglaplace_excess_kurtosis(c: float) -> float:
    c = _validate_c(c)
    if c <= 4:
        return np.nan
    c2 = c * c
    num = 6 * (
        2 * c**10 + 138 * c**8 - 615 * c**6 - 449 * c**4 - 132 * c2 - 24
    )
    den = c2 * (c2 - 16) * (c2 - 9) * (2 * c2 + 1) ** 2
    return num / den


def loglaplace_entropy(c: float) -> float:
    c = _validate_c(c)
    return 1.0 + np.log(2.0 / c)


for c in [0.8, 1.5, 2.5, 5.0]:
    print(
        f"c={c:>4}: mean={loglaplace_mean(c)}, var={loglaplace_variance(c)}, "
        f"skew={loglaplace_skewness(c)}, excess_kurt={loglaplace_excess_kurtosis(c)}"
    )

c= 0.8: mean=nan, var=nan, skew=nan, excess_kurt=nan
c= 1.5: mean=1.8, var=nan, skew=nan, excess_kurt=nan
c= 2.5: mean=1.1904761904761905, var=1.3605442176870748, skew=nan, excess_kurt=nan
c= 5.0: mean=1.0416666666666667, var=0.10540674603174603, skew=3.0046141326124665, excess_kurt=40.717785467128024

5) Parameter Interpretation

The shape parameter \(c\) is best understood via the log transform:

\[\log X \sim \mathrm{Laplace}(0, 1/c).\]

• Larger \(c\) means smaller Laplace scale \(1/c\) in log-space, so \(X\) concentrates more tightly around \(1\).

• Smaller \(c\) means heavier tails on the original scale (more extreme small and large values).

A quick rule of thumb from the tail:

\[\mathbb{P}(X > x) \approx \tfrac12\,x^{-c}\quad (x\ge 1).\]

So \(c\) acts like a Pareto tail index on the right tail.

6) Derivations

6.1 Expectation (and when it exists)

Using the piecewise PDF and splitting at 1:

\[\mathbb{E}[X^k] = \int_0^\infty x^k f(x;c)\,dx = \frac{c}{2}\int_0^1 x^{k+c-1}\,dx + \frac{c}{2}\int_1^\infty x^{k-c-1}\,dx.\]

The first integral is finite when \(k>-c\); the second is finite when \(k<c\). Evaluating gives:

\[\mathbb{E}[X^k] = \frac{c}{2}\Bigl(\frac{1}{k+c} + \frac{1}{c-k}\Bigr)=\frac{c^2}{c^2-k^2}, \quad -c<k<c.\]

Setting \(k=1\) shows the mean exists iff \(c>1\).

6.2 Variance

For \(c>2\) we have \(\mathbb{E}[X]\) and \(\mathbb{E}[X^2]\), so

\[\mathrm{Var}(X)=\mathbb{E}[X^2]-\mathbb{E}[X]^2\]

which simplifies to

\[\mathrm{Var}(X)=\frac{c^2(2c^2+1)}{(c^2-4)(c^2-1)^2}.\]

6.3 Likelihood (iid sample) and MLE

For data \(x_1,\dots,x_n\) with \(x_i>0\) in the standard model (only \(c\) unknown),

\[f(x_i;c)=\frac{c}{2x_i}\exp\bigl(-c|\log x_i|\bigr).\]

The log-likelihood is

\[\ell(c)= n\log c - n\log 2 - \sum_{i=1}^n \log x_i - c\sum_{i=1}^n |\log x_i|.\]

Differentiating and setting to zero:

\[\ell'(c)=\frac{n}{c}-\sum_{i=1}^n |\log x_i|=0\quad\Rightarrow\quad \hat c=\frac{n}{\sum_i |\log x_i|}.\]

Interpretation: since \(\log X\) is Laplace, \(|\log X|\) is Exponential(rate \(c\)), so inference on \(c\) reduces to inference on an exponential rate.

def loglaplace_loglik_c(x: np.ndarray, c: float) -> float:
    """Log-likelihood for the standard model (c only)."""
    x = np.asarray(x, dtype=float)
    c = _validate_c(c)
    if np.any(x <= 0):
        return -np.inf

    n = x.size
    s = np.sum(np.abs(np.log(x)))
    return n * np.log(c) - n * np.log(2.0) - np.sum(np.log(x)) - c * s


def loglaplace_mle_c(x: np.ndarray) -> float:
    x = np.asarray(x, dtype=float)
    if np.any(x <= 0):
        raise ValueError("all observations must be > 0")
    s = np.sum(np.abs(np.log(x)))
    if s <= 0:
        raise ValueError("degenerate sample: sum |log x| must be > 0")
    return x.size / s


# Quick check: simulate and verify the MLE is close to truth
c0 = 2.5
x = loglaplace_rvs_numpy(c0, size=10_000, rng=rng)
print("c0     =", c0)
print("c_hat  =", loglaplace_mle_c(x))

c0     = 2.5
c_hat  = 2.4890218861144686

7) Sampling & Simulation

7.1 Inverse CDF method (NumPy-only)

Draw \(U\sim\mathrm{Unif}(0,1)\) and transform using the quantile function:

\[\begin{split} X= \begin{cases} (2U)^{1/c}, & U<\tfrac12,\\ \left(\frac{1}{2(1-U)}\right)^{1/c}, & U\ge\tfrac12. \end{cases} \end{split}\]

This is exactly what loglaplace_rvs_numpy implements.

7.2 Alternative viewpoint (log transform)

Sample \(Z\sim\mathrm{Laplace}(0,1/c)\) and set \(X=\exp(Z)\). This is conceptually helpful, but the inverse-CDF above is simpler to implement without relying on a Laplace RNG.

8) Visualization

We’ll visualize:

• the PDF and CDF for different \(c\),

• Monte Carlo samples (including a log-scale view).

# Grids that resolve both sides around x=1
x_left = np.logspace(-3, 0, 400, endpoint=False)
x_right = np.logspace(0, 3, 400)
xgrid = np.unique(np.concatenate([x_left, x_right]))

c_values = [0.8, 1.5, 3.0, 7.0]

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))

for c in c_values:
    fig.add_trace(
        go.Scatter(x=xgrid, y=loglaplace_pdf_std(xgrid, c), mode="lines", name=f"c={c}"),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=xgrid,
            y=loglaplace_cdf_std(xgrid, c),
            mode="lines",
            name=f"c={c}",
            showlegend=False,
        ),
        row=1,
        col=2,
    )

fig.update_xaxes(title="x", type="log", row=1, col=1)
fig.update_xaxes(title="x", type="log", row=1, col=2)
fig.update_yaxes(title="density", row=1, col=1)
fig.update_yaxes(title="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380, legend_title_text="shape c")
fig.show()

# Monte Carlo: samples on x-scale and log-scale
c_mc = 2.5
n_mc = 50_000
x_samp = loglaplace_rvs_numpy(c_mc, size=n_mc, rng=rng)
z_samp = np.log(x_samp)

zgrid = np.linspace(np.quantile(z_samp, 0.001), np.quantile(z_samp, 0.999), 500)
laplace_pdf = 0.5 * c_mc * np.exp(-c_mc * np.abs(zgrid))  # Laplace(0, 1/c_mc)

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("Samples on original scale (log x-axis)", "log(samples) vs Laplace(0,1/c)"),
)

# Original-scale samples (histogram) with theoretical PDF overlay
fig.add_trace(
    go.Histogram(x=x_samp, nbinsx=120, histnorm="probability density", name="samples", opacity=0.6),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(x=xgrid, y=loglaplace_pdf_std(xgrid, c_mc), mode="lines", name="theory"),
    row=1,
    col=1,
)

# Log-scale samples should look Laplace
fig.add_trace(
    go.Histogram(x=z_samp, nbinsx=120, histnorm="probability density", name="log samples", opacity=0.6),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(x=zgrid, y=laplace_pdf, mode="lines", name="Laplace pdf"),
    row=1,
    col=2,
)

fig.update_xaxes(title="x", type="log", row=1, col=1)
fig.update_xaxes(title="z = log x", row=1, col=2)
fig.update_yaxes(title="density", row=1, col=1)
fig.update_yaxes(title="density", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()

# Monte Carlo check of mean/variance when they exist (c>2)
print("theory mean:", loglaplace_mean(c_mc))
print("MC mean    :", x_samp.mean())
print("theory var :", loglaplace_variance(c_mc))
print("MC var     :", x_samp.var())

theory mean: 1.1904761904761905
MC mean    : 1.1881560637980781
theory var : 1.3605442176870748
MC var     : 0.9558356968621015

9) SciPy Integration (scipy.stats.loglaplace)

SciPy provides a ready-to-use distribution object:

• pdf, logpdf, cdf, ppf

• rvs for sampling

• fit for MLE of parameters

We’ll use fit in two ways:

1. Shape-only fit with loc=0, scale=1 fixed (matches the “standard” model).

2. Full c, loc, scale fit (often less stable; interpret carefully).

from scipy.stats import loglaplace

c0 = 2.0
x = loglaplace.rvs(c0, size=5, random_state=rng)
print("SciPy rvs:", x)

# Evaluate
print("pdf:", loglaplace.pdf(x, c0))
print("cdf:", loglaplace.cdf(x, c0))

# Fit: generate data from the standard model
data = loglaplace_rvs_numpy(c0, size=5_000, rng=rng)

# Closed-form MLE for the standard model
c_hat_closed = loglaplace_mle_c(data)

# SciPy fit with loc/scale fixed to match the standard model
c_hat_scipy, loc_hat, scale_hat = stats.loglaplace.fit(data, floc=0.0, fscale=1.0)

print("c0          =", c0)
print("c_hat closed=", c_hat_closed)
print("c_hat SciPy =", c_hat_scipy)
print("(loc,scale) =", (loc_hat, scale_hat))

# Full fit (c, loc, scale) — can be sensitive for heavy tails
c_fit, loc_fit, scale_fit = stats.loglaplace.fit(data)
print("full fit (c,loc,scale)=", (c_fit, loc_fit, scale_fit))

SciPy rvs: [2.60896 0.72759 0.92112 0.64573 1.02715]
pdf: [0.05631 0.72759 0.92112 0.64573 0.92279]
cdf: [0.92654 0.26469 0.42424 0.20849 0.52608]
c0          = 2.0
c_hat closed= 2.004309521405735
c_hat SciPy = 2.004309521405735
(loc,scale) = (0.0, 1.0)
full fit (c,loc,scale)= (2.0137861206612757, -0.0034857193644913936, 1.000027126069282)

10) Statistical Use Cases

10.1 Hypothesis testing

In the standard model, inference on \(c\) is especially simple.

Since \(|\log X|\sim \mathrm{Exp}(\text{rate}=c)\), the statistic

\[S = \sum_{i=1}^n |\log x_i|\]

satisfies

\[2cS \sim \chi^2_{2n}.\]

This gives exact confidence intervals and tests for \(c\).

10.2 Bayesian modeling

The standard-model likelihood for \(c\) is proportional to

\[c^n\,\exp\bigl(-cS\bigr),\]

so a Gamma prior on \(c\) is conjugate:

\[c\sim\mathrm{Gamma}(a,b) \;\Rightarrow\; c\mid x \sim \mathrm{Gamma}(a+n,\; b+S).\]

(Here \(b\) is a rate parameter.)

10.3 Generative modeling

If you want multiplicative heavy-tailed noise, you can model

\[Y = \theta \cdot X, \quad X\sim\mathrm{LogLaplace}(c),\]

which is equivalent to adding Laplace noise in log-space:

\[\log Y = \log \theta + Z, \quad Z\sim\mathrm{Laplace}(0,1/c).\]

from scipy.stats import chi2

# Exact CI and an exact two-sided test for c in the standard model
c0 = 2.0
n = 300
x = loglaplace_rvs_numpy(c0, size=n, rng=rng)
S = np.sum(np.abs(np.log(x)))

alpha = 0.05
ci = (
    chi2.ppf(alpha / 2, df=2 * n) / (2 * S),
    chi2.ppf(1 - alpha / 2, df=2 * n) / (2 * S),
)

c_hat = loglaplace_mle_c(x)
test_stat = 2 * c0 * S  # under H0, this is chi2_{2n}
p_left = chi2.cdf(test_stat, df=2 * n)
p_two_sided = 2 * min(p_left, 1 - p_left)

print("c0   =", c0)
print("c_hat=", c_hat)
print("95% exact CI for c:", ci)
print("two-sided exact p-value for H0: c=c0:", p_two_sided)

c0   = 2.0
c_hat= 2.0630561150287168
95% exact CI for c: (1.836183726795856, 2.3029522418961106)
two-sided exact p-value for H0: c=c0: 0.6061428639658745

# Conjugate Bayesian update for c (standard model)
# Prior: c ~ Gamma(a0, rate=b0)
a0, b0 = 2.0, 1.0

c0 = 2.0
n = 200
x = loglaplace_rvs_numpy(c0, size=n, rng=rng)
S = np.sum(np.abs(np.log(x)))

a_post = a0 + n
b_post = b0 + S

posterior_mean = a_post / b_post
posterior_ci = (
    stats.gamma.ppf(0.025, a=a_post, scale=1 / b_post),
    stats.gamma.ppf(0.975, a=a_post, scale=1 / b_post),
)

print("prior (a,b)=", (a0, b0))
print("posterior (a,b)=", (a_post, b_post))
print("posterior mean=", posterior_mean)
print("posterior 95% CI=", posterior_ci)

prior (a,b)= (2.0, 1.0)
posterior (a,b)= (202.0, 93.21649390394049)
posterior mean= 2.166998473555129
posterior 95% CI= (1.8784506146804234, 2.4758606604180544)

11) Pitfalls

• Invalid parameters: \(c\le 0\) (and scale <= 0 in SciPy) is invalid.

• Nonexistent moments: mean requires \(c>1\), variance requires \(c>2\), etc. For smaller \(c\), Monte Carlo estimates like the sample mean can look “unstable” because the target quantity is infinite.

• Extreme samples: the inverse CDF produces very large values when \(U\) is extremely close to 1. In floating point code, clip uniforms away from 0/1.

• Numerical evaluation: prefer logpdf for likelihoods to avoid underflow, especially when \(x\) is tiny/huge or \(c\) is large.

• Fitting with loc/scale: unconstrained fit may be sensitive for heavy-tailed data; consider fixing loc/scale when you have a clear standardization.

12) Summary

• Log-Laplace models positive, heavy-tailed data and is Laplace in log-space.

• The PDF is piecewise power-law and the right tail behaves like a Pareto tail with index \(c\).

• Raw moments satisfy \(\mathbb{E}[X^k]=\frac{c^2}{c^2-k^2}\) for \(|k|<c\), so some moments may not exist.

• Sampling is easy with an inverse CDF (NumPy-only).

• In the standard model, \(|\log X|\) is exponential, giving a closed-form MLE, exact chi-square inference, and a conjugate Gamma posterior for \(c\).